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We quantify the consistency of numerical-relativity black-hole-binary waveforms for use in 
gravitational-wave (GW) searches with current and planned ground-based detectors. We compare 
previously published results for the {I = 2, |m| = 2) mode of the gravitational waves from an equal- 
mass nonspinning binary, calculated by five numerical codes. We focus on the lOOOM (about six 
orbits, or 12 GW cycles) before the peak of the GW amplitude and the subsequent ringdown. We 
find that the phase and amplitude agree within each code's uncertainty estimates. The mismatch 
between the {I — 2, \m\ = 2) modes is better than 10"'' for binary masses above &QMq with respect 
to the Enhanced LIGO detector noise curve, and for masses above 180 Mq with respect to Advanced 
LIGO, Virgo and Advanced Virgo. Between the waveforms with the best agreement, the mismatch 
is below 2 x 10~*. We find that the waveforms would be indistinguishable in all ground-based de- 
tectors (and for the masses we consider) if detected with a signal-to-noise ratio of less than ~ 14, 
or less than ~ 25 in the best cases. 

PACS numbers: 



I. INTRODUCTION 

Direct detection of gravitational waves is expected in 
the next few years by a network of ground-based laser- 
interferometric detectors, LIGO [J, [1,11], Virgo [1, H] and 
GEO which operate in the frequency range ^ 

10^-10^ Hz. The scientific scope of gravitational-wave ob- 
servations will be extended (see Q for a recent overview) 
by space-based instruments such as LISA [H, [ll[ , which 
will be sensitive to signals at significantly lower frequen- 
cies. A likely source for the first detection, and an es- 
sential part of the science objectives of all gravitational- 
wave detectors, is the merger of black-hole-binary sys- 
tems. Detection of gravitational-wave events and their 
further analysis rely on the theoretical modeling of wave- 
forms. Until recently, theoretical waveforms for the co- 
alescence of black holes were based on analytic approx- 
imations to the full general theory of relativity, in par- 
ticular the post-Newtonian expansion, which models the 
signal from the slow inspiral of the two black holes, and 
black-hole perturbation theory, where the complex ring- 
down frequencies of black holes can be computed (see 



[T3j for reviews). These methods cannot currently 
model from first principles the merger phase, when the 
wave amplitude peaks. Correspondingly, data analysis 
methods so far had to be developed without information 
from complete black-hole-binary waveforms. 

The situation changed with breakthroughs in numer- 
ical relativity in 2005 [H, [H, [11] that made it possible 
to calculate the late inspiral, merger and ringdown of a 
black-hole-binary system in full general relativity, and to 
calculate the gravitational waves produced in the pro- 
cess. Since that time many more numerical simulations 
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and efforts have been made 
to produce waveform templates based on numerical re- 
sults m, [tI, [tI, [11] . Some of these template banks are 
already available to be used for searches within the LSC 
Algorithm Library [sij . There is also an ongoing project 
to test search pipelines with injections of numerical data 
into simulated LIGO and Virgo noise, the Numerical IN- 
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Jection Analysis (NINJA) project [s^]- The work de- 
scribed in this paper was conceived as complementary 
to that in the NINJA effort, and has subsequently been 
dubbed the Samurai project. 

If waveform template banks based (at least in part) 
on numerical results are to be confidently used in detec- 
tor searches, it is important to know the accuracy of the 
input numerical waveforms. Most numerical waveforms 
are published with some internal error analysis and un- 
certainty estimates. However, our goal in this paper is 
to perform a stronger consistency check by comparing 
the results of different numerical codes, produced using 
different formulations of the Einstein equations, initial 
data, gauge conditions, and numerical techniques. First 
studies comparing numerical waveforms were performed 
in [H, [131; our project extends that earlier work. 

This paper serves two purposes: (1) to verify that 
the numerical waveforms we compare agree with each 
other within the uncertainty estimates originally pub- 
lished with those waveforms, and (2) to quantify the 
differences between the waveforms in terms and mea- 
sures meaningful to both the numerical-relativity and 
gravitational- wave data-analysis communities. In par- 
ticular, in addition to making a direct comparison be- 
tween the phase and amplitude of the respective wave- 
forms, we also compute their mismatch with respect to 
the Enhanced LIGO, Advanced LIGO, Virgo and Ad- 
vanced Virgo detectors, and the maximum signal-to-noise 
ratio (SNR) below which these waveforms would be in- 
distinguishable . 

In the present comparison wc focus on the physi- 
cal system that has been studied in most detail by 
the numerical-relativity community: a binary consisting 
of equal mass, non-spinning black holes following non- 
eccentric inspiral and merger. Specifically, we consider 
the dominant (i = 2, |m| = 2) spherical-harmonic mode, 
which has been the focus of most data-analysis research 
to date, and is the most important from the point of view 
of GW detection. It is now known that the sub-dominant 
modes are also important for parameter estimation, but 
we will leave an anlysis of those to future work; see also 
our comments in the Conclusion. In order to keep the 
data analysis aspects of this paper straightforward, we 
will discuss only single detectors, and neglect the sub- 
tleties introduced when dealing with networks of detec- 
tors or with time-delay interferometry. 

We will consider the waveform from roughly lOOOA/ be- 
fore merger, to about 80A/ after merger, where M is the 
total binary mass in geometrical units. This is about six 
orbits before merger, or 0.005 {M/Mq) seconds, where 
Mq — 1.477 X 10'^ m is the mass of the sun. 

We make use of results from the BAM [13, Isl. CCATIE 
[Slj . Hahndol [H, [83] and MayaKranc codes |37|. which 
all use the BSSN/moving-puncture [H, [H, [s^,!!, HI, S 
HHjHII approach and finite-difference techniques, and the 
SpEC code , which solves a variant of the generalized- 
harmonic system [s^, [H, [131 using pseudospectral meth- 
ods. 



The paper is organized as follows. In Section |TT] we 
summarize the numerical waveforms that we analyze, and 
the codes that were used to produce them. In Section Hill 
we directly compare the phase and amplitude of the wave- 
forms. In Section HVl we calculate the detector mismatch 
between the waveforms for a range of masses and detec- 
tors, and determine the SNR below which the waveforms 
would be indistinguishable for a single GW detector. In 
Section |V] we draw some conclusions from our compar- 
isons. 



II. NUMERICAL WAVEFORMS 

A. The physical system 

Wc restrict our attention to the modeling of one phys- 
ical system, the orbital inspiral of two black holes of 
equal mass and zero spin with vanishing eccentricity. 
The size of the orbits decreases as the system loses en- 
ergy through gravitational radiation emission, until the 
black holes merge to form a single spinning black hole. 
The final mass and spin, which determine for exam- 
ple the ringdown frequencies, are routinely determined 
from numerical simulations with a variety of methods, 
for early results see e.g. [13, [H, [H, [H, [o^. Cur- 
rent simulations are able to obtain very accurate re- 
sults for the final black hole parameters, and as exam- 
ples we quote the results and error estimates reported 
for the SpEC and BAM codes [s^ [6l| , which are consistent 
within the given error estimates: the mass of the final 
black hole has been found as M/ = 0.95162 ± 0.0000271/ 
with SpEC, and as Mf = 0.9514 ± 0.0016M with BAM; 
the dimensionless spin (Kerr parameter) has been com- 
puted as Sf/M'^ = 0.68646 ± 0.00004 with SpEC and as 
Sf/Mf = 0.687 ± 0.002 with BAM. This corresponds to 
a dominant ringdown frequency of the £ = m = 2 mode 
of Mu = 0.5539 ± 0.0018 (conservative BAM estimate) or 
Mw = 0.5535 ± 0.00003 (SpEC estimate), using interpo- 
lation in tabulated values of ringdown frequencies given 
in [93. 

If we decompose the gravitational-wave signal from 
this system into spherical harmonics, the (£ = 2,m — ±2) 
modes dominate. The frequency of these two modes 
(which are related by a tt/2 phase shift) is very close 
to twice that of the orbital motion during inspiral, and 
steadily increases as the black holes approach merger. 
The amplitude of the signal is a function of the fre- 
quency, and also increases. The signal frequency and 
amplitude peak around merger, and then the ampli- 
tude decays exponentially as the merged black hole rings 
down to a stationary Kerr black hole. For an equal- 
mass nonspinning binary, we know from numerical simu- 
lations that the peak frequency is approximately fpeak ~ 
16(Mq/M) kHz. Five orbits before merger, the wave fre- 
quency is / w 1.95(MQ/Af) kHz, and one hundred orbits 
before merger the frequency is / w O.38(M0/M) kHz, as 
estimated by post-Newtonian methods. 
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FIG. 1: Theoretical noise curves (power spectral density Sn, 
see Eq. I12p for the detectors Enhanced LIGO, Advanced 
LIGO, Virgo and Advanced Virgo. 
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FIG. 2: The gravitational-wave strain from an optimally- 
oriented 60 Mq equal-mass nonspinning black-hole binary lo- 
cated 100 Mpc away from the detector. The waveform covers 
about six orbits, or twelve GW cycles, before merger. 



The most sensitive ground-based detectors currently 
in operation, LIGO and Virgo, can detect signals from 
black-hole binary (BH-BH) mergers out to distances of 
up to several hundred Mpc (dependin g on the binary's 
mass and orientation, see e.g., [l|, [tqI. l98l|V Estimated 
event rates for BH-BH coalescence events are on the order 
of one e very few year s, but with large uncertainties , see 
e.g. l9§,[lM[lpil,ll0l. The future Enhanced LIGO jlO^ 
and Virgo-I- 1] detectors will increas e eve nt rates by ~ 5 
times, whe reas the Advanced LIGO |104| and Advanced 
Virgo [l05l | detectors will increase event rates by roughly 
three orders of magnitude as compared with the current 
detectors. 

These detectors are sensitive to frequencies ranging 
from ^ 10 — 40 Hz up to ~ 2 kHz. The merger signal will 
be in this frequency range for systems with total masses 
of roughly 5-250 Mq. The merger will be in the most sen- 
sitive part of the detectors' frequency bands for masses 
around 50 M©, and for that case the detectors will also 
be sensitive to the signal from the last ten orbits before 
merger. Theoretical estimates of the noise curves for the 
four detectors we consider, Enhanced LIGO, Adva nced 
LIGO Virgo [H and Advanced VIRGO [lo|, are 
shown in Figure [T] (for Virgo and Advanced LIGO we use 
approximate analytical formulas as displayed in [79j). 

Astrophysical black hol es rnay form b inaries through 
a number of mechanisms [lOOl . Il06l . Il07l | . In general the 
black holes will have different masses and wi ll be spin- 
ning (high spins may be typical [I08l . ll09l . IllOj ). and the 
orbits will be eccentric. But gravitational-radiation emis- 
sion reduces the eccentricity, so for typical comparable- 
mass inspirals the eccentricity is expected to be negligible 
by the time the binaries have reached the frequen- 
cies we have just dis cusse d (for the situation in globular 
clusters see however jl07l |). For this reason most analyt- 
ical and numerical work in modelling gravitational-wave 
signals has focussed on binaries that follow non-eccentric 
(or "quasi-circular" ) inspiral; see however [H, 113, [H, [fill 
for numerical results on eccentric binaries. 



The preceding discussion motivates a focus on the last 
orbits before merger of binaries following non-eccentric 
inspiral. We study a binary that consists of black holes 
with equal mass and no spin, simply because this config- 
uration has been studied in the most detail in numerical- 
relativity simulations. We consider the gravitational- 
wave signal from the last ~ 6 orbits and merger of 
this system. Figure [2] shows one polarization of the 
gravitational-wave strain from an example of such a bi- 
nary, with total mass 60 Mq, optimally oriented to the 
detector and located 100 Mpc away. 



B. Numerical codes 



To calculate the gravitational-wave signal in full gen- 
eral relativity, we first require a solution of Einstein's 
equations. This must be produced numerically, and can 
be done in a number of different ways. We will compare 
the results of five computer codes: BAM, CCATIE, Hahndol, 
MayaKranc and SpEC. These differ in their procedure for 
constructing black-hole-binary initial data, decomposi- 
tion of Einstein's equations into a numerically well-posed 
and stable form, numerical techniques used to evolve the 
data, choice of gauge conditions during the evolution, and 
details of the calculation of the gravitational-wave signal. 

We will summarize the different methods of setting up 
the initial data, the formulations of Einstein's equations, 
and the numerical techniques. The purpose is not to pro- 
vide a full exposition of these methods (the full technical 
details can be found in the references given in Table |l| 
but to make clear the similarities and differences of the 
five codes. 
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1. Initial data 

Four of the codes, BAM, CCATIE, MayaKranc and 
Hcihndol, use Bowen-York puncture initial data jll2| . 
The chief features of these data are that the spatial 
metric is conformally flat, and the physical momenta 
and spins of the black holes can be specified directly as 
parameters in t he B owen-York solution of the momen- 
tum constraint |ll3l |. The black holes manifest them- 
selves in the data through topological wormholcs, which 
also allow the spatial slices to bypass the black-hole 
singularities. The wormholes are compactified so that 
their ends are mapped to single points, or "punctures" 
[iil[ni[ni[iil,[n3- The nature of these wormholes 
changes during a dynamical evolution [9l|, IllSl . IllOl . Il20l | , 
but the spatial slices never reach the singularities, and 
the data can be evolved without recourse to excising any 
region of the computational domain. 

The data used in each code differ only in the choices of 
initial separation of the black-hole punctures, and their 
momenta, and in the method of numerically solving the 
Hamiltonian constraint in the puncture approach; this 
last distinction will affect the accuracy of the solution of 
the Hamiltonian constraint, but we assume in this work 
that the solution used by all four codes is sufficiently 
accurate that the remaining numerical errors do not con- 
tribute to the differences we measure in the final results. 

Ideally the momenta are chosen to produce non- 
eccentric quasi-circular inspiral. In the simulation from 
the Halmdol code, the initial momenta were modified by 
hand until a roughly quasi-circular inspiral was obtained. 
The BAM, CCATIE, and MayaKranc simulations used pa- 
rameters calculated by post-Newtonian methods, as out- 
lined in [56j . This procedure results in an eccentricity of 
e < 0.0016. The choices of initial momenta are given in 
Table HI 

The SpEC code uses excision data: the data extend only 
to the black-hole apparent horizons. The data were con- 
structed by sol ving the c onformal-thin-sandwich initial- 
value equations [l2ll Il22l | , with suitable boundary condi- 
tions on the apparent horizons and at the outer bou ndary 
to produce non-spinning black holes [l23l . Il24l . Il25j | in an 
orbit with small radial velocity [3^ . The parameters ap- 
propriate to quasi-circular inspiral were first predicted by 
the methods described in [l25j . and then modified using 
the iterative procedure described in [3l|, H^j , to yield an 
eccentricity of below e ^ 5 x 10~^. Once again the data 
are conformally flat. 

The reader should not make too much (or too 
little) of the differences between these two types 
of data, Bowcn- York-puncture and conformal-thin- 
sandwich-excision. Although they are constructed in 
quite different ways, both sets of data are based on sim- 
ilar choices of the free data in initial-value equations (in 
particular conformal flatness), and may not be physically 
very different; elucidating their exact differences is not 
trivial. Conversely, they are not identical, and there is 
no reason to expect a priori that the waveforms resulting 



from evolutions of both sets of data will precisely agree. 
Evaluating that difference is part of the analysis in this 
work. 



2. Evolution systems 

The codes that start with pimcture initial data — BAM, 
CCATIE, MayaKranc and Hahndol — evolve the data with 
the BSSN formulation of Einstein's equations, and follow 
[isl . [l6t in the use of coordinate conditions that allow the 
black holes to move across the grid. The BSSN evolution 
system [sl, [s^ [qO] is combined with hyperbolic evolu tion 
equations for the lapse and shift (1-l-log slicing [l26l | and 
the f-driver shift condition [13, [sO]), which have been 
shown to lead to a well-posed initial- value proble m |93 | . 

The SpEC code uses a first-order formulat ion IMI of 
the generalized- harmonic-gauge system |93, Il27| with 
built-in constraint-damping terms |l28l . Il29| | . This sys- 
tem is manifestly symmetric hyperbolic and well-posed. 
The gauge source functions are chosen to be constant 
in a comoving frame during the inspiral p^ . and are 
evolved according to a sourced wave equation during 
merger and ringdown [6l|. The characteristic fields of 
the system are all outward-flowing (into the holes) at 
the excision boundaries, so no boundary conditions are 
nee ded or imp osed there. The outer boundary conditions 
jol . Il30l . Il3l| are desi gned to p revent the influx of con- 
straint violations [T3^ . Tl33l . [iM [Tssl . [Hi, ITstI. [l38l and 
undesired incoming gravitational radiation [l3ij | . while 
allowing the outgoing gravitational radiation to pass 
through the boundary. 



3. Numerical techniques 

The moving-puncture codes solve the partial- 
differential equations of the BSSN formulation of the 
Einstein equations with finite-difference methods. The 
numerical domain consists of nested Cartesian domains, 
such that successive levels of refinement are placed both 
around the individual black holes (and centered on 
the punctures) and around the entire black-hole-binary 
system (centered on the origin of coordinates). The 
details of the mesh refinement differ between codes, and 
the full details can be found in the relevant references. 
The spatial finite-differencing is in general fourth-order 
in Hahndol, CCATIE, and MayaKranc, and sixth-order 
in BAM. Here, centered differences are used with the 
exception of shift advcction terms, which use one-point 
lopsided stencils. Integration forward in time is per- 
formed with a fourth-order Runge-Kutta method. The 
Hcihndol code uses a uniform time step, while the other 
BSSN codes use variants of a Berger-Oliger scheme, 
where finer grids can evolve with smaller time steps. 
The refinement boxes that are centered around the black 
holes move with them through the grid. 
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The SpEC code uses multidomain pseudospectral meth- 
ods on a grid with two excised regions, one just inside 
the apparent horizon of each hole, and employs a dual- 
frame technique to track the motion of the holes [l^ . The 
computational domain consists of two sets of concentric 
spherical shells, one surrounding each excised region, an- 
other set of concentric spherical shells extending to the 
outer boundary, and a structure of touching cylinders 
that fills in the remaining volume and overlaps some of 
the spherical shells. Inter-domain bo unda r y co nditions 
are enforced with a penalty method |l40l . Il4l| . Time 
stepping is accomplished via the method of lines, using 
an adaptive fourth/fifth order Runge-Kutta method. 



4- Summary of the numerical codes and waveforms 

Table U summarizes the similarities and differences of 
the five waveforms, and the codes used to produce them. 
More details on the waveforms can be found in the fol- 
lowing references: The BAM waveform is from the highest- 
resolution D12 simulation described in The CCATIE 
results have been obtained from the simulation described 
in [4^. The MayaKranc simulation is the e = simula- 
tion described in [s^ . The SpEC waveform corresponds to 
the waveform "30c-l/N6" described in Ref. [6l|; the in- 
spiral portion of this waveform is more comprehensively 
discussed in Ref. [3^, including a detailed error anal- 
ysis. The Hahndol waveform comes from the highest- 
resolution (grid spacing of M/32 at the finest level) evo- 
lution of the "di = 10.8 M" data presented in [H]. For 
all simulations the nominal Courant factor was 0.5, al- 
though for the BAM, CCATIE and MayaKranc simulations 
the Courant factor was lowered on the two outermost 
mesh-refinement levels. 

No new simulations were performed for this paper, al- 
though the MayaKranc waveform results from an updated 
extrapolation procedure as described in Section [II CI 

Table [J also provides uncertainty estimates in the GW 
phase and amplitude, quoted separately for the inspiral 
regime (up to a frequency of roughly Mu; = 0.2), and the 
merger and ringdown regime. It is important to bear in 
mind that each uncertainty estimate applies only to the 
waveform, and not to the code used to produce it. For ex- 
ample, a code that uses second-order-accurate finite dif- 
fercing may well produce waveforms more accurate than 
any presented here, if run at sufficiently high resolution 
with sufficiently accurate initial data, and if the gravita- 
tional waveforms were extracted sufficiently far from the 
source. 

Note that the apparent accuracy of the phase and am- 
plitude depend strongly one how one chooses to align 
waveforms from different simulations, and whether quan- 
tities are considered as functions of time, phase or fre- 
quency. All of these choices are valid when comparing 
results produced by evolving the same initial data with 
the same evolution system, and varying only numerical 
resolution, radiation extraction radii and outer bound- 



ary location, and one is free to make the choice that 
gives the lowest error estimate. As such, the methods 
used to estimate the phase and amplitude errors differ 
for each waveform; more details can be found in some of 
the references listed in Table ID 

Having said that, in the present study we are compar- 
ing results from different codes, with different initial data 
and gauge conditions, and the disagreements we see from 
different waveform alignment choices may exaggerate, or 
hide, the "real" differences between the waveforms. 

We find that the least ambiguous method of compari- 
son is to plot quantities with respect to the frequency Muj 
of the {£ = 2,m = 2) mode of ^'4. This choice removes 
the need to apply a time and phase shift when compar- 
ing the wave amplitude, and the freedom of a constant 
phase shift in a phase comparison is straightforward to 
interpret. 



C. Extraction of gravitational waves 

In the numerical simulations presented here, the grav- 
itational waves are extracted using the N ewman-Penrose 
Weyl tensor component ^4 jl42l ! Il43| . which at infi- 
nite separation from the source is related to the complex 
strain h = /i+ — i/ix by 14J|, 



h = lim / dt' [ dt"^'4. 

'■^°°Jo Jo 



(1) 



Note that the amplitude of the gravitational-wave 
strain falls off as 1/r, where r is the distance of the de- 
tector (or, in a numerical code, the extraction sphere) to 
the source, and so we generally consider rh (and r^!^), 
which in the weak-field region will be independent of r. 

It is useful to discuss gravitational radiation fields in 
terms of spherical harmonics of spin- weight s = —2, 
'^in-n which represent symmetric tracefree 2-tensors on 
a sphere, and in this paper we will only consider the 
dominant i = 2, m = ±2 modes, with basis functions 



Y^-:\ = xl--{l-cosere-''^, 



2 „-2i< 



647r 



-2 _ 



-(l + cos.)^e- 



(2) 



i.e., we will consider the cases £ = 2,m ^ ±2 of the 
projections 

h,,„ = {Yf^ ,h)= r ThY^ smOdO, (3) 
Jo Jo 

of the complex strain h (bar denotes complex conjuga- 
tion). In the nonspinning case considered here, we have 
equatorial symmetry so that h22 = h2-2, and 

h(i) = yj^e^'* ((1 + C0s9f h22(i) + (1 - COS^)' h22(t)^ 
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TABLE I: Summary of numerical codes. The initial separation is the coordinate separation between the punctures (for moving- 
puncture codes) or between the centers of the excision surfaces (SpEC). The initial momenta specified in the moving-puncture 
codes are (pt.Pr) /M, where pt/M is the tangential momentum and Pr/M is the radial momentum. The SpEC parameters are 
described in [32|. "Bulk FD order" indicates the spatial finite difference order in the bulk of the computational domain (i.e., 
not including mesh-refinement boundary zones), /imin is the spatial resolution on the finest mesh-refinement level or domain. 
The wave extraction radii are given, and r^x —> oo indicates that the results were then extrapolated to infinity. The references 
provide full details of the implementation of the codes and the simulations that were used in this study. For the CCATIE result no 
numerical convergence results were published, but based on the code specification and resolution for this run, a phase accuracy 
between the Hahndol and BAM/MayaKranc results can be assumed. The amplitude errors quoted for CCATIE and MayaKranc were 
estimated for the present paper and were previously unpublished. 
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(xlO-^*) 
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uncertainty 
(radians) 


uncertainty 
(percentage) 














insp. 


merger 


insp. 
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Fiiutc-difference moving-puncture codes 


BAM [84, 85] 


D = 12M; (0.085, -5.373 x 10"'') 


6 


19 


90M 


e < 0.0016 


0.1 


1.0 


4.0 


6.0 


CCATIE [511 


D = IIM; (0.090, -7.094 x 10"'') 


4 


20 


120Af 


e < 0.0016 






2.0 


5.0 


Hahndol [25, 86, 87] 


D = 10.8M; (0.0912,0.0) 


4 


19 


60A/ 


e < 0.008 


2.4 


5.0 


10.0 


10.0 


MayaKranc [37] 


D = 12M; (0.085, -5.343 x 10"'') 


4 


15.5 


rex oo 


e < 0.0016 


0.1 


1.1 


4.0 


8.0 


Pscudospectral excision code 


SpEC [31, 61J 


D = 14.436Af; 


n/a 


~ 3 


rex oo 


e < 5 x 10"^ 


0.006 


0.02 


0.1 


0.3 


rc^^ = 0.41360Af, MQq = 0.016708, 
Vr^-4.26 X 10"", /r = 0.939561 



















The coordinate radius at which the waves were ex- 
tracted from the numerical solution is given for each code 
in TablelH For the MayaKranc and SpEC codes, the waves 
were extracted at several radii, and then extrapolated to 
r^x — > cxD, to give a more accurate estimate of the wave 
that would be measured by a distant GW detector. The 
extrapolation procedure involves aligning the waveforms 
with respect to some definition of retarded time [3l| , and 
then treating the error due to extraction radius as a poly- 
nomial in powers of 1 /r^x (Sll . [sst . Different polynomial 
fits were performed for the inspiral and merger for the 
MayaKrainc waveform, and the specific extrapolation pro- 
cedure used for the SpEC waveform is given in [sil, [6l| . 
Waves were extracted from the CCAT IE s imulation using 
the Zcrilli-Moncricf procedure (see [l45j | for a review), 
from which ^'4 can be readily derived. 

The direct waveform comparisons in Section IIIII deal 
with . The data-analysis comparisons in Section IIVI 
are based on the strain, rh. To produce the strain from 
^4 one needs "merely" to integrate twice with respect 
to time, as in Eq. ([1]), and choose appropriate integra- 
tion constants. However, this procedure is not as trivial 
as it at first appears. One might naively assume that 
integration constants could be chosen on simple physi- 
cal grounds, for example that the strain rings down to 
zero after the black holes have merged, and that it oscil- 
late around zero at all times. Such requirements have 
been found to work adequately in some cases for the 
{£ = 2,m = 2) mode, but even in the best cases unusual 
artifacts remain, and these become more pronounced 
when one considers higher modes; see [146| for some ex- 
amples. One reason for these difficulties is that the wave- 



forms contain small numerical errors and gauge effects, 
which become greatly exaggerated when integrated over 
the entire duration of the waveform — and to calculate 
the strain we must perform such an integration twice. 

A tempting alternative is to work only in the Fourier 
domain. Start with the numerically generated ^'4(t), cal- 
culate the Fourier transform, ^'4(7), and then it is trivial 
to perform two time integrations to obtain the Fourier 
transform of the strain, 



h(/) 



^4(/) 

47r2/2- 



(4) 



The integration constants have been ignored in this pro- 
cedure, or, rather, they have been implicitly set to zero. 
If we now perform an inverse Fourier transform to cal- 
culate h{t), we will recover similar artifacts to those we 
would have seen if we had performed two time integra- 
tions of *4(t). 

However, in this paper we use this very method to 
calculate h(/) to use in our match calculations. Our jus- 
tification is that for a selection of waveforms we have 
independently calculated h(t) by a number of different 
methods (with varying levels of success in removing nu- 
merical and gauge artifacts in the final strain), and have 
then used the Fourier transform of this quantity in match 
calculations, and obtained very similar results; we will 
discuss the impact of the small differences that we see 
in Section IIVI Our conclusion is that the choice of in- 
tegration constants, and modifications that "clean" the 
waveform of non-physical artifacts, although they may 
lead to serious differences in the time-domain waveform, 
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do not significantly afi'ect the match calculation for the 
{e = 2,m = 2) mode. 



III. DIRECT COMPARISON OF PHASE AND 
AMPLITUDE 

We now compare the waveforms produced by the five 
codes. For the purposes of gravitational-wave detection, 
the most meaningful comparison will include the noise 
spectrum of the detector. We will make comparisons rele- 
vant to detection and parameter estimation in Sect ion ITVl 
In this section we directly compare the numerical wave- 
forms in a manner that is independent of any particular 
detector. The quantities we will compare are the ampli- 
tude A{t) and the phase (j){t) of the {£ = 2,m = 2) mode 
of r\E'4, which are defined by 

r*4,22(i) = A(t)e-'*(*). (5) 

The GW frequency for the (£ = 2, m = 2) mode is given 
byc^(<) 

The amplitude and phase are the two pieces of raw 
output from the computer code that define the wave- 
form, so they allow the most direct comparison between 
results from different codes. More generally, an ampli- 
tude/phase comparison allows us to quantify waveform 
differences independent of any detector — if two wave- 
forms accumulate one cycle of phase disagreement during 
the last ten cycles before they reach the peak amplitude, 
that is a difference that will exist no matter which detec- 
tor they pass through. 

On the other hand, there are a number of ambigui- 
ties in an amplitude/phase comparison, which we will 
describe as we proceed. For the purpose of gravitational- 
wave detection, the detector mismatch and SNR are more 
meaningful quantities to compare. We can summarize 
the situation as follows: a direct comparison of ampli- 
tude and phase is most useful to the numerical relativist, 
while the mismatch and SNR are most useful to the data 
analyst. 



A. Phase 

In comparing the phases of two waveforms from dif- 
ferent simulations, (pii^i) ^tnd 4>2it2), we cannot simply 
calculate (/>i(ii) — ^2(^2), because the time coordinates ti 
and t2 may not be the same. The two simulations may 
have been started at different points along the binary in- 
spiral, meaning that ti = does not label the same event 
as <2 = 0. More simply: although the waveforms may be 
identical, one will reach a detector later than the other. 

To make a comparison we first have to decide on an 
event at which the two waveforms should agree, and to 
then apply a relative time shift and phase shift so that 
the chosen event occurs at the same time and phase for 
each waveform. For example, if we were to align the 
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FIG. 3: The GW frequency as calculated from the raw BAM 
data, and as given by the fitting procedure described in the 
text. The two lines are indistinguishable if viewed in black- 
and-white. 



phases at the time when the waveform amplitude reaches 
a maximum, then we would first determine the times 
Ti and T2 when each waveform's amplitude reaches its 
maximum, and then study the quantity 



A(/.(t) = <t)i[t + Ti)-<t,2{t + T2) + <i>[T2) - c^{Ti), (6) 



where by construction the amplitude maxima now occur 
at t = and A0(O) = 0. 

The problem with this procedure is that A(/)(t) is ex- 
tremely sensitive to the accuracy with which Ti and T2 
were determined, particularly around the merger, when 
the GW frequency increases rapidly. One solution is to 
make a further small time shift, until the overall phase 
disagreement between the two waveforms has been min- 
imized. Such a suggestion has been used in the past in 
matching NR and FN waveforms [78l.[l47t. and in NR-FN 
comparisons [g^- This, however, is an approach designed 
not to determine the differences between two waveforms, 
but to minimize them. Another option, which avoids the 
time-shift ambiguity altogether, is to compare the phases 
as a function of GW frequency (pi^Mij)] this procedure 
was used in [25j . and we will use it here. 

The GW frequency as read from the numerical data is 
too noisy at early and late times to allow a clean direct 
parametric plotting of phase vs frequency. We instead fit 
the frequency t o a c ombination of the TaylorT3 FN fre- 
quency formula |l48l | and a modification of the frequency 
ansatz introduced in Specifically, the TaylorTS ex- 
pression for the orbital frequency of the binary during 
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inspiral is given up to 3.5PN order by |l48l . Il49| 
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where we have given the last (3.5PN) term an arbitrary 
coefficient, a{i'), where v is the symmetric mass ratio 
ly = TOim2/M^. This term is known in PN theory, but 
we will instead fit it to our numerical data, given our 
modified definition of the variable r, which we will now 
discuss. The definition of t in standard PN theory is 



iy{tc - t) 
5M '■ 



(8) 



where tc is a PN estimate of the "coalescence time" . The 
expression ([7]) diverges when r = 0, so in order to pro- 
duce a formula which can be fit through our data, we use 
instead 



25Af2 



1, 



(9) 



an d we now treat tc as a parameter to fit to the data, as 
in [l48j |. With our new definition of r, the expression ([7]) 
becomes inaccurate near r = (which is anyway true for 
any post-Newtonian expression near merger), but does 
not diverge. To model the ringdown phase, we modify 
the ansatz suggested in [g^ , and write the full frequency 



n(t) = npNiT) + 

{Qf-npN{T)) 



1 + tanh[ln ~ [t - to)/b]Y 



(10) 



The constants {tc, to, k, a, b, flj} are parameters that are 
determined to produce the best fit to the numerical data. 
The constant flf corresponds to a fit of the ringdown 
frequency. The frequency as a function of time is shown 
for the BAM code in Fig. [31 as calculated from the raw 
numerical data, and as given by the fitting procedure we 
have just described; the GW frequency is related to the 
orbital frequency by a factor of two. The GW phase as 
a function of frequency for each of the five waveforms is 
shown in Fig. 21 




FIG. 4: GW phase as a function of frequency Muj, for 
the five codes. The frequency is given both in terms of the 
dimensionless orbital frequency Muj, and the frequency in kHz 
scaled with respect to the total mass of the binary in solar 



Figure |5l compares the phase of each waveform with 
that from the SpEC code. The GW frequency lOOOM 
before merger, where our waveforms nominally begin, 
is close to Mui = 0.055, and this is the frequency at 
which our comparison begins. After merger, the merged 
black hole rings down to the Kerr solution, and the GWs 
are emitted at the ringdown frequency, which is close 
to Mlo = 0.55; this is where we end our comparison. 
(The precise ringdown frequency for the equal-mass, non- 
spinning, zero-eccentricity configuration is Mlo ~ 0.5535 

ill-) 

In the left panel of Figure we show the phase dis- 
agreement during inspiral, ending at Mlo = 0.2, which is 
reached about half an orbit before merger. A phase shift 
is applied so that the phases all agree at Muj = 0.055. 
We sec that the accumulated phase disagreement is be- 
low 0.3 radians for all codes except Hahndol, for which 
the larger eccentricity and lower resolution lead to larger 
dephasing against the SpEC results. The behaviour of 
the three other waveforms is roughly consistent with the 
numerical methods used to produce them: the BAM wave- 
form was produced with the highest-order spatial finite- 
differencing (sixth-order), and while fourth-order spatial 
finitc-diffcrencing was used to produce both the CCATIE 
and MayaKranc results, the MayaKranc simulation was 
performed at slightly higher resolution, and the results 
were further extrapolated with respect to radiation ex- 
traction radius. 

The most important point is that the results of each 
code agree within their respective uncertainty estimates. 

The right panel of Figure |5l shows the accumulated 
phase disagreement during the last orbit, merger and 
ringdown. The phases are shifted to agree at the low- 
est frequency shown in the figure, Muj = 0.2, so that 
we can see how the phase disagreement behaves during 
the merger regime only. Note that the waveform from 
the Hahndol simulation becomes very noisy late in the 
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FIG. 5: Phase comparison. The left panel shows the phase comparison between the SpEC waveform and the others during 
inspiral, from Mlj = 0.055 up to Mtj = 0.2, which is about one orbit before merger. The corresponding uncertainty estimates are 
0.1 rad (BAM and MayaKranc) and 2.4 rad (Hahndol). The right panel shows the phase comparison during merger and ringdown, 
from A/cij = 0.2 up to A/a; = 0.55. The uncertainty estimates during merger are, in radians, 1.0 (BAM), 1.1 (MayaKranc), and 5.0 
(Hahndol). In both panels a phase shift was applied so that the phases agreed at the lowest frequency shown. 



ringdown, which accounts for the poor behaviour above 
AIu! « 0.52. Note also that while the merger and ring- 
down plot sweeps through roughly twice the range of 
frequencies as the inspiral plot, the length of time cov- 
ered during the inspiral (about 900M) is much greater 
than that during the merger (about 180M). In this 
sense the phase disagreement grows more quickly during 
merger. The phase disagreements of the different wave- 
forms are again consistent with uncertainty estimates, 
but are larger for merger and ringdown than they were 
during inspiral. 



The two panels of Fig. [S] were designed to show sep- 
arately the phase difference accumulated during inspi- 
ral. or during mergcr/ring-down. When considering the 
phase as a function of frequency (as done in Fig. [5|), the 
only freedom is an overall additive constant to the phase. 
Thus, the total accumulated phase difference during in- 
spiral and merger/ringdown can be obtained by vertically 
offsetting the curves in the right panel of Fig. O so that 
the phase-differences at Muj = 0.2 agree in both panels. 
For instance, the total accumulated phase-difference be- 
tween BAM and SpEC at Mlu = 0.52 would be the sum 
of 0.1 rad (from the left panel of Fig. O, and 0.28 rad 
(from the right panel), i.e., 0.38 rad. For the other 
codes, one finds at Muj = 0.52 the following total ac- 
cumulated phase-differences relative to SpEC: Hahndol 
1.36 rad, CCATIE 0.55 rad, and MayaKranc 0.18 rad. 
The reader may choose to calculate the total accumu- 
lated phase disagreement at any frequency, although one 
should bear in mind that beyond Mlo = 0.52 the curves 
in Fig. O are less reliable, due to errors in the curve fit 
Eqn (|10|) through noisy numerical data. 



B. Amplitude 

In comparing the GW amplitude between codes, we 
once again consider the amplitude as a function of fre- 
quency, A{Muj), which is shown for the five codes in 
Fig. [6l The amplitude comparison is shown in Fig. [T] 
Once again the comparison during inspiral, shown in the 
left panel, covers the frequency range Mlu G [0.055,0.2], 
and the comparison during merger and ringdown, shown 
in the right panel, covers the frequency range Mlo G 
[0.2,0.55]. The oscillations in the right panel of Fig. [7] 
are probably due to small errors from gauge effects that 
are exaggerated in this plot by the rapid change in the 
GW frequency near merger. Recall that the Hahndol 
waveform becomes unreliable at about Afo; 0.52. 

Note once again that the agreement is within the esti- 
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FIG. 6: The amplitude as a function of GW frequency, 
A{uj) — |r*I'4,22| for the five codes. 
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mated uncertainties of the waveforms. 

The conclusion of our direct comparison of the GW 
phase and amphtude is that all five codes are consistent 
within their stated error bars. This provides an impor- 
tant consistency check on the numerical accuracy and 
validity of each waveform. Not only that, it provides 
us with an upper limit on the variations in the wave- 
forms due to different choices of initial data (in these 
cases puncture data versus quasi-equilibrium conformal- 
thin-sandwich excision data), and different gauge choices. 
The latter can lead to noticeable differences in the ampli- 
tude and phase of the waveform extracted from the simu- 
lation (see, for example, the discussions in [6ll. [sl . [l50t ) . 
With suitable gauge choices these differences should de- 
crease as the waves are extracted successively further 
from the source, and indeed it is usually possible to per- 
form some procedure to extrapolate the waveform to es- 
timate the result that would be measured infinitely far 
away [Sll, [H, [6ll, [6^. These are delicate procedures, 
and one may still worry that the different gauge choices 
between codes will lead to large differences in the final 
waveforms. In this section we have shown that, if such 
differences exist, they are small and within the error bars 
of each simulation. 

The results so far provide information that allow nu- 
merical relativists to quantify the accuracy and consis- 
tency of their results. In the next section we will make 
comparisons relevant to data analysis and GW astron- 
omy. 



IV. DETECTION 

A more meaningful comparison from the point of view 
of GW detection is th e best match (and mismatch) be- 
tween waveforms [l5l| . 

The match is usually calculated in the frequency do- 
main. Consider two time series x{t) and y(t), which will 
be the two waveforms wc wish to compare. The Fourier 
transform is given by 



Snif) m 



x{t)e^-'''f' dt . 



(11) 



(In LSC applications the opposite sign convention is used 
for the phase in the Fourier transform definition, but the 
choice of sign does not affect the results here.) In prac- 
tice we calculate a discrete Fourier transform on the nu- 
merical data. We calculate the time when the wave am- 
plitude reaches its maximum, tmax, and then truncate 
the waveform lOOOM before this time, and 80Af after. 
The resulting truncated waveform is then resampled ev- 
ery O.IM, to give a data set with 10,800 points. We then 
take a discrete Fourier transform of each such data set, 
and retain only the half of the data set that covers posi- 
tive frequencies. We also verified that our results did not 
change significantly when the sampling rate was varied. 

We can define an inner product between i(/) and 
y{f) weighted with the noise spectrum of the detector. 



{x\y) := 4 Re 



Snif) 



df 



(12) 



In the same way we define a norm of a waveform x(f) by 
\x\ = \/{x\x). 

The signal-to-noise ratio (SNR) is defined with respect 
to this waveform norm. Recall that throughout this pa- 
per we have been dealing with rh and r\['4, where r is the 
distance of the detector (or numerical wave extraction) 
from the source, and one should remember to use the real 
strain h in the definition of the SNR. For clarity, let us 
define h = rh, as calculated from the numerical code, and 
then the SNR is given by 



R 



(13) 



where R is the distance of the source from the detector, 
usually in units of Mpc. 

The best match [l5l| is defined as the inner product 
{x\y) normalized by the norms of each waveform, and 
maximized over relative time and phase shifts (r and <&) 
between the two waveforms: 



M = max ^ '^ -^ 



(14) 



We can view this procedure as adjusting the waveforms 
with respect to their time of arrival, and their initial 
phase, such that we achieve the best agreement. 

One convenient way [l53| | to calculate the best match 
with respect to the phase shift is to first normalize 
each polarization of the two waveforms as ei+_x = 
i+,x/|a;+,x| and 62+, x = y+^x /|y+,x |, and to define 

A = (ei+|e2+)2 + (el+|e2x)^ 
B = (elx|e2+)^ + (elx|e2x)^ 
C = (ei+|e2+)(eix|e2+) + (ei+|e2x)(eix|e2x)- 

In general one should also orthornomalize the two wave- 
forms, but in this work wc consider only optimally ori- 
ented binaries, and so this is not necessary. The best 
match is then given by [l53j | 



= max 



A + B 



A-B 



(15) 



For the waveforms we consider here the match is very 
close to unity, and it makes more sense to quote the mis- 
match, defined as 1 — M.. 

In evaluating p^ . we must choose /min and /max- Ide- 
ally these would be (0, 00), but in practice they are based 
on the range of frequencies for which the Fourier trans- 
form is reliable (very high and very low frequencies con- 
tain unphysical artifacts due to numerical errors and the 
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FIG. 7: Comparison of the amplitude as a function of GW frequency, A{u}). The left panel shows the percentage disagreement 
during inspiral (up to A4uj — 0.2). The corresponding amplitude uncertainties are 2% (CCATIE), 4% (BAM and MayaKranc), and 
10% (Hahndol). The right panel shows the same quantity during merger and ringdown, for which the uncertainties are 5% 
(CCATIE), 6% (BAM), 8% (MayaKranc), and 10% (Hahndol). The vertical dashed line indicates the approximate location of the 
amplitude maximum, Muj = 0.5. 
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FIG. 8: The mismatch between the SpEC waveform and each of the other codes (the line colors match those in previous plots). 
The three plots show the results for the Enhanced LIGO, Advanced LIGO, Virgo and Advanced Virgo noise curves. The lower 
end of the mass range was chosen such that the entire numerical waveform was included in the detector's frequency band. 



sampling rate of the data), and the relevant frequency 
window of a given detector. The frequency range for the 
Enhanced LIGO detector is chosen as 30 Hz up to 2 kHz, 
and for the Advanced LIGO, Virgo and Advanced Virgo 
detectors it is from 10 Hz up to 2 kHz. The acceptable 
frequency range for the Fourier transforms of the numer- 
ical data is from /M = 0.001 up to /M = 1. The actual 



integrals are performed over the intersection of the de- 
tector and waveform frequency ranges. 

The time shift r deserves some discussion. From our 
experience with phase comparisons in Section IHII we 
know that in general we must time-shift two waveforms 
with respect to each other in order to realistically esti- 
mate their agreement. It is natural to apply that time 
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shift to one of the waveforms in the time domain, but 
if we have the entire frequency-domain representation of 
the signal, we can also calculate the effect of a time shift 
on the match in the frequency domain. In our case we do 
not have the entire frequency-domain waveform: the nu- 
merical waveform was truncated lOOOM before merger, 
and 80M after merger, and as such the calculated Fourier 
power is incorrect outside a certain range of frequencies. 
Put another way, a time shift in the time domain would 
result in using in our analysis a different 1080M-long por- 
tion of one of the waveforms, and there is no way that 
the analysis of the Fourier transform of one 1080A/-long 
segment of the waveform can capture the effect of choos- 
ing a different lOBOM-long segment. As such we apply a 
time shift to one of the waveforms before calculating its 
Fourier transform, and maximise the match with respect 
to this time shift. 

Figure [8] shows the minimum mismatch between the 
SpEC waveform and each of the other four, for the En- 
hanced LIGO, Advanced LIGO, Virgo and Advanced 
Virgo noise curves. A mass range of 60 — 300 M© was 
used for the Enhanced LIGO detector, and 180- 300 Mq 
for the others. The lower mass was dictated by the desire 
that the waveform begin below the low-frequency cut- 
off of the detector. At yet lower masses we would need 
to use longer waveforms, i.e., waveforms that extend to 
lower frequencies. We can see from these results that the 
mismatches are excellent: below 10"'^ in all cases, and 
for all except the Hahndol waveform, the mismatch is 
below 4 X lO^''. The performance of each waveform is 
consistent with the expected accuracy of each code, and 
with the results presented in Section Hill Recall also that 
the mismatch calculation usually involves adjusting the 
mass associated with one of the waveforms, in order to 
improve the result. If such a minimisation (or, in terms 
of the match, a maximisation) were performed here, the 
mismatches would improve further. 

The mismatch is in general very sensitive to phase dif- 
ferences, and this may explain the worse mismatch of 
the Hahndol waveform, which shows large variations in 
the phase disagreement, due mostly to higher eccentric- 
ity and lower numerical resolution (see Fig. [5l where the 
larger dephasing due to is apparent). We should note 
however that the mismatch is still extremely small, and 
easily meets the standard detection criteria, which we 
will discuss below. 

As we pointed out in Section Hi C[ the Fourier trans- 
form of the strain used for match calculations was pro- 
duced from the Fourier transform of ^^4, i.e., we calcu- 
lated 5*4 (t) ^4(7) — ^ h(/). The matches we calculate 
are very close to unity, and we wish to know how much 
the results vary if we first calculate the strain h{t) in 
the time domain, or if we vary slightly the length of the 
waveforms in either time or frequency. We tested the 
robustness of some of our match calculations to these 
changes, and found that the results could vary by as much 
as 2 X 10~^, but the values shown in Fig. [5] were almost 
always lower than those calculated by other methods. As 



such, we consider the curves in Fig. [8] as lower bounds on 
the mismatch, and note that we expect that in the worst 
case they would be no more than 2 x 10~^ higher. 

The best mismat ch r e quired for detection is usually 
given as 0.035 [H, [HI |lli, [l56l. f or which no more 
than 10% of signals will be lost 15l| . This is the best 
mismatch required between a member of a waveform tem- 
plate bank and the true physical waveform that is de- 
tected. In LIGO detector searches, templates are con- 
structed such that the worst mismatch betwe en succes - 
sivc members of the template bank is 0.03 [98l . ll54l . Il55l |. 
This places the more stringent requirement on the accu- 
racy of the theoretical waveform of a mismatch better 
than 0.005; see the discussion in [iSTj . This threshold is 
well above the largest mismatches calculated here. The 
conclusion, then, is that current numerical waveforms are 
sufficiently accurate for detection purposes with all cur- 
rent and planned ground-based detectors. 



A. Parameter estimation 

We now evaluate the differences between the wave- 
forms with respect to measurement of the source's pa- 
rameters. The theory of para meter estimatio n acc uracy 
is developed and discussed in jl52l . Il58l . Il59l . Il60l | . We 
defer a detailed analysis of these waveforms with respect 
to parameter estimation to future work, and here mak e 
a first analysis based on the criterion proposed in |157| . 

Imagine that a GW signal is detected, and the wave- 
forms studied here are used to estimate the parameters of 
the source. Each of our numerical waveforms is slightly 
different, and if the detected signal were strong enough 
(i.e., the SNR were large enough) the estimated parame- 
ters of the source would be different depending on which 
numerical waveform we use. However, if the SNR is be- 
low a certain value, any two of our waveforms will be 
indistinguishable. To put this discussion in context, a re- 
liable detection requires an SNR above a threshold which 
is usually chosen between about 5 and 8, depending on 
details of the det ector and search performed (compare 
e.g. m, [Tm, [HI), Therefore, if the SNR has to be be- 
low this threshold value for two waveforms to be indis- 
tinguishable, then it is meaningless to claim that they 
agree: we will never be able to perform an experiment to 
check. 

We can estimate the highest SNR for which two wave- 
forms are indistinguishable for a single GW detector as 
follows. Choose the binary mass, the detector, and dis- 
tance of the source to the detector. Determine the time 
and phase shift such that the mismatch between the two 
waveforms is a minimum. Calculate the difference be- 
tween those two aligned waveforms in the time domain, 
6h{t) = hi(i) — h2{t), and transform to the frequency do- 
main to produce 5h{f). It was shown in [l57j | that when 
the inner product of ^h(/) satisfies the criteria 



(dh\Sh) < 1, 



(16) 
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the two waveforms are indistinguishable. The left-hand 
side of Eqn pB]) and the SNR are both inversely pro- 
portional to the distance of the source to the detector. 
Therefore, having calculated the SNR and the value of 
((5h|(5h) for one source distance, we may immediately es- 
timate the maximum SNR such that the inequality in 
Eqn. ([TBI) is satisfied. The results of this calculation arc 
shown in Figure [HI 

This analysis applies to any estimation of the intrinsic 
parameters of the binary, like the total mass and mass 
ratio. But we should emphasize that, since we minimized 
((5h|(5h) with respect to a phase and time shift, Fig.[9]does 
not apply to an estimation of the signal's time of arrival 
or phase, or parameters that rely on them (like the sky 
location). It should also be emphasized that our analysis 
is restricted to parameter estimation based on the output 
of only one detector. 

Note also that the relative performance of each wave- 
form is not necessarily the same for both the detection 
and measurement analyses. This is because the match 
calculation finds the best agreement in phase only, while 
the measurement calculation locates the best match in 
phase and amplitude. 

We can estimate from Figure IH] that the SpEC, BAM, 
CCATIE and MayaKranc waveforms will be indistinguish- 
able according to Eqn. [16] if the SNR is below about 25; 
if the Hahndol waveform is to also be indistinguishable, 
the SNR must be below 14. 

SNRs of above 25 are expected to be uncommon for 
the Enhanced LIGO and Virgo detectors; for example, 
in the NINJA study numerical waveforms were injected 
into simulated detector noises at SNRs no higher than 
30 [S^l- For the Advanced LIGO and Virgo detectors, 
however, which have roughly ten times the range, SNRs 
in excess of 25 arc far more likely. 



V. CONCLUSION 

We have compared numerical-relativity waveforms for 
the last six orbits, merger and ringdown of an equal- 
mass nonspinning binary with minimal eccentricity, as 
produced by five different computer codes. We focussed 
on the (I = 2,m = 2) mode. The purpose was to per- 
form a stringent consistency check of the results from 
these codes. We verified that accuracy in the waveform 
phase and amplitude for each code was consistent with 
the uncertainty estimates originally published with each 
waveform. 

We also calculated the best mismatch between the 
most accurate waveform (calculated with the SpEC code) 
and each of the others, for the Enhanced LIGO, Ad- 
vanced LIGO, Virgo and Advanced Virgo detectors, and 
found that it was below 10~^ in all cases, and below 
2 X lO"** in the best cases. Recall that the criteria such 
that no more than 10% of signals is lost is 0.005 (assum- 
ing a standard template-bank spacing). 

Finally, we calculated the maximum SNR below which 



the signals would be indistinguishable if observed in the 
Enhanced LIGO, Advanced LIGO, Virgo or Advanced 
Virgo detectors. For the best cases this is about 25, and 
is never lower than 14. This suggests that these numer- 
ical waveforms are more than adequate as ingredients in 
template banks for GW searches and parameter estima- 
tion (of intrinsic parameters, at least) with the Enhanced 
LIGO and Virgo detectors, for which an SNR above 25 is 
unlikely. The Haindol waveform is distinguishable from 
the others at an SNR of only 14, in which case a more de- 
tailed study would be necessary to compare its parameter 
estimation fidelity with the other waveforms. Nonethe- 
less, we estimate that less accurate waveforms would not 
be desirable for GW data analysis purposes. We expect 
that these results extend to other numerical waveforms, 
if they exhibit similar or better levels of numerical un- 
certainty. 

An important caveat to the above statements is that 
these waveforms could only be used "as is" in detec- 
tor searches for high-mass binaries. Detection of bina- 
ries with lower masses would require waveforms that are 
longer (extend to lower frequencies), for example by com- 
bining analytic approximations (usually post-Newtonian 
and effective-one-body [EOB] waveforms) with numer- 
ical results. Methods have been pro pose d to produce 
both hybrid waveforms [13, [zl, H l80l.ll47l| and analytic 
waveforms based on either a phenomenological ansatz 
[zl, III, HSl or the adjustment of free parameters in var- 
ious EOB prescriptions [H, H [H, [11, [TH, [HI. It is 
the accuracy of those "complete" waveforms that will 
be important in lower-mass searches, and one may also 
find that sufficiently accurate complete waveforms will re- 
quire either much longer numerical waveforms as input, 
or more physically accurate approximation techniques. 
Such questions are beyond the scope of this paper, but 
are an important topic for future work. 

Furthermore, for general waveforms higher spherical 
harmonic modes will play a much more important role, in 
particular for parameter estimation. T his h as first been 
pointed out for ground-based detectors [l65j using post- 
Newtonian inspiral waveforms, with much recent work 
on ground-based detectors Il66ll a s well as on the pla.nned 
space-based LISA mission [10llll|. see e.g. 
Recently significant improvements in parameter estima- 
tion for LISA from higher mode contri butions ha ve also 
been hinted at for numerical waveforms |l70l . Il7l| . Large 
values of the SNR will be typical for future generations of 
ground based detectors, and even more so for LISA detec- 
tions. This will make it possible to determine source pa- 
rameters far more accurately than with current ground- 
based detectors, which will in turn place more stringent 
accuracy requirements on numerical waveforms. How- 
ever, we hope that by the time LISA flies (2018-I-), and 
by the time second and third generation ground based 
interferometers are in operation, the typical accuracy of 
numerical waveforms will have far surpassed that of those 
considered in this study. Our more immediate concern 
is whether current numerical codes are producing wave- 
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FIG. 9: The signal-to-noise ratio (SNR) below which the SpEC and each other waveform will be indistinguishable in any 
measurement of parameters. Results are shown for the Enhanced LIGO, Advanced LIGO, Virgo and Advanced Virgo detectors. 
See text for further explanation. 



forms of sufficient numerical and physical accuracy for 
use in current data-analysis applications, and our results 
suggest that they are, at least as far as one is concerned 
with the quadrupole mode, which is typically the basis 
of current matched-filter searches. 
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